# Fig 4c. - Map of SBS40b in Czech Republic counties

rm(list = ls())

# Spatial objects
library(geojsonio)
library(sp)
library(sf)
library(raster)
# Data handling
library(dplyr)
library(stringr)
# Plots
library(ggplot2)
library(ggspatial)
library(ggsflabel)
library(ggnewscale)
library(ggpattern)

# Function for plotting histogram
plot_histo <- function(df,var,namevar=var,viridis=F,alpha=1){
  df_plot <- filter(df, get(var) != "Missing")
  p <- ggplot(df_plot, aes_string(x = var, fill = var)) +
    geom_bar(stat = "count", color = "black", width = 1,alpha=alpha) +
    stat_count(geom="text", colour="white", size=3.5,
               aes(label=..count..), position=position_stack(vjust=0.5)) +
    xlab(namevar) + ylab("Number of samples") +
    theme_classic() +
    theme(legend.position = "none", axis.text.x = element_text(angle=45, vjust=1, hjust=1))
  if (viridis){p <- p + scale_fill_viridis_d()}
  p
}

# Loading individual data
geocode_longest_CZ <- readRDS("datasets/geocode_longest_CZ")
summary_CZ_all <- readRDS("datasets/summary_CZ_all")

# Loading maps
# World countries
world_countries <- geojson_read("maps/stanford-kk522dt9425-geojson.json", what = "sp")
# Subsetting European countries
europe_countries <- world_countries %>% subset(region_un == "Europe")

# Regions from Czech Republic
cz_region <- readRDS("datasets/cz_region")

# CZ borders
cz_country=aggregate(cz_region)

# CZ regions where centers are located
cz_centers <- cz_region %>% subset(ref %in% c("JC", "JM", "OL", "ST", "PR"))
cz_centers@data <- cz_centers@data %>% transmute(ref2 = ifelse(ref == "PR", "ST", ref))
cz_centers <- cz_centers %>% aggregate(by="ref2")

# CZ minicharts information
cz_region_minichart <- cz_region %>% subset(ref %in% c("JC", "JM", "OL", "PR"))
cz_region_minichart@data <- cz_region_minichart@data %>% mutate(center_region = case_when(Region_Code == "PHA" ~ "Prague region",
                                                                                          Region_Code == "JHM" ~ "Brno region",
                                                                                          Region_Code == "JHC" ~ "Ceske Budejovice \n region",
                                                                                          Region_Code == "OLK" ~ "Olomouc region"))

# Theme for plot
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_rect(fill = "aliceblue"),
                         plot.background = element_rect(fill="white"),
                         panel.border = element_blank(),
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         axis.title.y = element_blank(),
                         plot.title = element_text(size=22)))

map.test <- ggplot() +
  # European countries contours & names
  geom_sf(data=st_as_sf(europe_countries), alpha=1, lwd=0.5, fill="ghostwhite")+
  # CZ border
  geom_sf(data=st_as_sf(cz_country), alpha=1, lwd=3)+
  # CZ regional RCC incidence
  geom_sf(data=st_as_sf(cz_region), aes(fill = Incidence_1999_2013), alpha=1, lwd=0,color=NA)+
  scale_fill_gradient("Kidney cancer incidence (ASR-W) \n (1998-2013)", low="#c6dbef",high="#08519c", 
                      guide = guide_colorbar(title.position ="left", barwidth=10, order=1)) + 
  # CZ borders of regions where the recruitment centers are localized
  geom_sf(data=st_as_sf(cz_centers), alpha=0, lwd=1)+
  new_scale_fill() +
  geom_sf_label_repel(data=st_as_sf(cz_region_minichart),aes(label=center_region), nudge_x=1.2,nudge_y=.3, size = 4.5)+
  # Individual points located inside CZ
  geom_spatial_point(data = filter(geocode_longest_CZ,county != "Slovakia"), crs = 4326, aes(x = jitter(lon_lgst,10), y = jitter(lat_lgst,10), fill = SBS1536A_cosmic_q), shape=21,size=5,alpha = 0.8) +
  scale_fill_viridis_d("SBS40b attribution \n (quartiles)", direction=-1, guide = guide_legend(title.position ="left")) +
  # Scale and north arrow
  annotation_scale(location = "bl", width_hint = 0.2, text_cex = 1) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  coord_sf(crs=4326,xlim = c(12.05, 18.9), ylim = c(48.5,51.1), expand = FALSE) +
  theme(legend.position = c(1,1), legend.justification = c("right", "top"), legend.direction = "horizontal", 
        legend.box.just = "right", legend.background = element_rect(fill="ghostwhite"), legend.key.height = unit(.9, "cm"),
        legend.text = element_text(size=12), legend.title = element_text(size=12)) +
  theme_opts

# Only keeping counties with strictly more than 10 cases 
df_barplot_CZ <- geocode_longest_CZ %>% 
  group_by(county2) %>% 
  filter(n() > 10) %>% ungroup()

bar.testplot_list <- 
  lapply(1:length(unique(df_barplot_CZ$county2)), function(i) { 
    gt_plot <- ggplotGrob(
      plot_histo(df_barplot_CZ[df_barplot_CZ$county2 == unique(df_barplot_CZ$county2)[i],], "SBS1536A_cosmic_q",alpha=0.5)+
        labs(x = NULL, y = NULL) + 
        scale_fill_viridis_d(direction=-1) +
        theme(legend.position = "none", rect = element_blank(),
              line = element_blank(), text = element_blank()) 
    )
    panel_coords <- gt_plot$layout[gt_plot$layout$name == "panel",]
    gt_plot[panel_coords$t:panel_coords$b, panel_coords$l:panel_coords$r]
  })
# Creating barcharts for each regions
bar_annotation_list <- lapply(1:length(unique(df_barplot_CZ$county2)), function(i) 
  annotation_custom(bar.testplot_list[[i]], 
                    xmin = summary_CZ_all$lon_modif[summary_CZ_all$county2 == as.character(unique(df_barplot_CZ$county2)[i])] - 0.2,
                    xmax = summary_CZ_all$lon_modif[summary_CZ_all$county2 == as.character(unique(df_barplot_CZ$county2)[i])] + 0.4,
                    ymin = summary_CZ_all$lat_modif[summary_CZ_all$county2 == as.character(unique(df_barplot_CZ$county2)[i])] - 0.1,
                    ymax = summary_CZ_all$lat_modif[summary_CZ_all$county2 == as.character(unique(df_barplot_CZ$county2)[i])] + 0.4) )

# Adding the barcharts to the main plot
result_plot <- Reduce(`+`, bar_annotation_list, map.test)
result_plot
# Output in 8x13
